function [psi_val,psi_x,psi_y] = boundary_function(x,y,Re,caseNum)
%function [psi_val,psi_x,psi_y] = boundary_function(x,y,Re,caseNum)
switch caseNum
    case 0  % for stream function with analytic solution
        [psi_val,psi_x,psi_y] = psi_h(x,y,Re,caseNum);
    case 1  % for square cavity
        psi_val = zeros(size(x));
        psi_x = psi_val;
        psi_y = psi_val;    
        
        idx = find((y-1.0)<1e-5);  % only the upper boundary's flux in y direction is 1
        psi_y(idx) = -1.0;
end

end